!Résolution du système d'équations différentielles de Jarek.
!Avec routine de tri à la fin (répartition homoigène de points dans gfx log-log).

module fonctions
implicit none
real, parameter :: y10 = 1. !condition initiale pour x[t]
real, parameter :: y20 = 1. !condition initiale pour y[t]
real, parameter :: gamma = 0.2 !rapport des vitesses gamma = c1/c2
!real, parameter :: gamma = 0.0 !rapport des vitesses gamma = c1/c2
real, parameter :: ti = 0. !temps initial
real, parameter :: tf = 11000. !gfx log-log
!real, parameter :: tf = 10.1 !gfx lin-lin
!real, parameter :: tf = 4.2 !3eme gfx
logical, parameter :: logscale = .true. !gfx log-log
!logical, parameter :: logscale = .false. !gfx lin-lin et 3eme gfx
real, parameter :: precis = 0.000001 !précision du résultat
integer, parameter :: itmax = 15 !nombre d'itérations maximal pour la convergence
integer, parameter :: nbeq = 2
contains

function fct(x,y)
real, intent(in) :: x
real, intent(in), dimension(:) :: y
real, dimension(size(y)) :: fct
	fct(1) = -4.*gamma*y(1)**2 -(3+gamma**2)*y(1)*y(2)
	fct(2) = -4.*y(2)**2 -(3+gamma**2)*y(1)*y(2)
return
end function fct
  
!------------------------------------------------------------------------
!Int‚gration de l'‚quation diff‚rentielle ordinaire:
! Dyk/Dx = f(x,y1,...,yN) ,yk(xi) = yik   ,k=1,N (C.I.) i=i=cte
!On doit donner:
!xi,xf : extr‚mit‚s de l'intervale d'int‚gration [xi,xf]
!yi : conditions initiales t.q. yik=y(xi). Tableau de dimension N
!itmax : nombre maximum d'it‚rations permis (a la fin de chaque it‚ration,
!        le nombre d'intervales est doubl‚
!eps   : erreur relative: tableau de dimension N avec eps(i) l'erreur
!        relative sur la i-Šme ‚quation diff‚rentielle
!
!On fait appel …:
!function fct(x,y)      : fonction qui pour chaque ensemble des valeurs des
!			  variables x et y = (y1,...,yN) fournit le tableau
!                         fct contenant les valeurs des fonctions
!			  f1(x,y),...,fN(x,y).
!
!Le programme sort les valeurs:
!yf    : valeur cherch‚e de y(xf): tableau contenant les valeurs des N
!	 variables d‚pendantes calcul‚es pour x=xf
!n     : nombre d'intervales qui a permis d'atteindre la pr‚cision d‚sir‚e
!x     : tableau de dimension (n+1) contenant les (n+1) valeurs interm‚diaires
!        de la variable x
!y     : tableau de dimension (N,n+1) contenant les valeurs des N ‚quations
!        diff‚rentielles aux n+1 points d'‚valuation.
!        on a un tableau de dimension (1:n+1) avec les valeurs 0,1,...,n
!
!ATTENTION: x,y sont … d‚clarer comme suit:
!             real, pointer, dimension(:) :: x
!             real, pointer, dimension(:,:) :: y
!           Le allocate est fait ici, mais il ne faut pas oublier de faire
!           deallocate(x,y) dans le programme principal.
!           Il faut toujours se rappeler que la d‚claration d'interface est
!           en bijection avec la d‚claration du programme principal

subroutine rk4sys(fct,xi,yi,xf,itmax,eps,yf,n,x,y)
integer, intent(out) :: n
integer, intent(in) :: itmax
real, intent(in) :: xi,xf
real, intent(in), dimension(:) :: yi,eps
real, intent(out), dimension(:) :: yf
real, pointer, dimension(:) :: x
real, pointer, dimension(:,:) :: y
interface
  function fct(x,y)
    real, intent(in) :: x
    real, intent(in), dimension(:) :: y
    real, dimension(size(y)) :: fct
  end function fct
end interface
integer, parameter :: ninit = 110000 !gfx log-log
!integer, parameter :: ninit = 200 !gfx lin-lin
!integer, parameter :: ninit = 50 !3eme gfx
integer :: npas, it, i
logical, dimension(size(yi)) :: masque
real, dimension(:), allocatable :: xt
real, dimension(:,:), allocatable :: yt
real :: h,xc
real, dimension(size(yi)) :: err,yc,yc_old,t1,t2,t3,t4

write (*,*) "Iteration         Points "
write (*,*) "-------------------------"
npas = ninit
it = 0
boucle: do
  it = it + 1
  if (it > itmax) then
    write (*,*) "rk4sys: pas de convergence"
    stop
  end if
  npas = 2*npas
  h = (xf-xi)/npas
  xc = xi
  yc = yi
  allocate (xt(npas),yt(size(yi),npas))
  do i=1,npas
    t1=h*fct(xc,yc)
    t2=h*fct(xc+h/2.,yc+t1/2.)
    t3=h*fct(xc+h/2.,yc+t2/2.)
    t4=h*fct(xc+h,yc+t3)
    yc = yc + (t1 + 2*t2 + 2*t3 + t4)/6.
    xc = xi + i*h
    xt(i) = xc
    yt(:,i) = yc
  end do
  if (it >= 2) then
    masque = .false.
!    err = abs((yc-yc_old)/(15.*yc))          !richardson relative
    err = abs((yc-yc_old)/15.)               !richardson absolue
    where (err <= eps) masque = .true.
    write (*,'("   ",i3,"       ",i10)') it,npas
    if (all(masque)) then
      yf = (16*yc-yc_old)/15.
      n = npas
      allocate(x(n+1),y(size(yi),n+1))
      x(1:n+1) = (/xi,xt(1:n)/)
      y(1:size(yi),1) = yi                    !conditions initiales
      y(1:size(yi),2:n+1) = yt                !solutions cherch‚es
      deallocate(xt,yt)
      exit
    end if
  end if
  yc_old = yc
  deallocate(xt,yt)
end do boucle
return
end subroutine rk4sys

end module fonctions

!----------------------------------------------------------------------
!----------------------------------------------------------------------

program ann2sp
use fonctions
implicit none
real, pointer, dimension(:) :: x
real, pointer, dimension(:,:) :: y
real, dimension(nbeq) :: eps,yi,yf
real :: xi,xf,h
integer :: n,i,un,j,nouveau

eps = precis
xi = ti
xf = tf
yi(1) = y10
yi(2) = y20
call rk4sys(fct,xi,yi,xf,itmax,eps,yf,n,x,y)
write(*,*)
write(*,*) "Saving the file..."
open(unit=1, file='ann2spx.txt')
open(unit=2, file='ann2spy.txt')
if (logscale .eqv. .false.) then
	i=1
	ecriture0: do
		if (x(i) < 1.0) then
			write(1,*) x(i),y(1,i)
			write(2,*) x(i),y(2,i)
			i = i+1
			if (i > n+1) exit ecriture0
		end if
		if ( (x(i) >= 1.0) .and. (x(i) < 2.)) then
			write(1,*) x(i),y(1,i)
			write(2,*) x(i),y(2,i)
			i = i+2
			if (i > n+1) exit ecriture0
		end if
		if ( (x(i) >= 2.) .and. (x(i) < 3.)) then
			write(1,*) x(i),y(1,i)
			write(2,*) x(i),y(2,i)
			i = i+3
			if (i > n+1) exit ecriture0
		end if
		if (x(i) >= 3.) then
			write(1,*) x(i),y(1,i)
			write(2,*) x(i),y(2,i)
			i = i+5
			if (i > n+1) exit ecriture0
		end if
	end do ecriture0
else !répartition logarithmique des valeurs
	cherche1: do i=2,n+1
		if ((x(i) >= 0.1) .and. (x(i-1) < 0.1)) then
			un = i
			exit cherche1
		end if
	end do cherche1
	h = log10(x(un+1)) - log10(x(un))
	write(1,*) x(un-1),y(1,un-1)
	write(1,*) x(un),y(1,un)
	write(2,*) x(un-1),y(2,un-1)
	write(2,*) x(un),y(2,un)
	j = un
	ecriture: do
		cherchei: do i=j+1,n+1
			if ( (log10(x(i)) >= h + log10(x(j))) .and. ( log10(x(i-1)) < h + log10(x(j))) ) then
				nouveau = i
				exit cherchei
			end if
			if (i==n+1) exit ecriture
		end do cherchei
		j = nouveau
		if (j > n+1) exit ecriture
		write(1,*) x(j),y(1,j)
		write(2,*) x(j),y(2,j)
	end do ecriture
end if
close(1)
close(2)
deallocate(x,y)
write(*,*)
write (*,*) "Program finished."
read(*,*)
end program ann2sp
